load("data/Wvs_analysis.Rdata")

###################################################################
###################3        TREATMENT     ##########################
###################################################################


wvs$natharbdist_km <- wvs$two_port_all_pred1/1000

wvs$natharbdist_km <- wvs$natharbdist_km/100

###################################################################
###################3        OUTCOMES     ##########################
###################################################################

#democracy is good thing
table(wvs$Q238)
wvs$democracy_good[wvs$Q238==1] <- 4
wvs$democracy_good[wvs$Q238==2] <- 3
wvs$democracy_good[wvs$Q238==3] <- 2
wvs$democracy_good[wvs$Q238==4] <- 1

#army rule is bad
table(wvs$Q237)
wvs$armyrule_bad <- wvs$Q237

#strong leader bad
wvs$strongleader_bad <- wvs$Q235

#theocracy bad
wvs$theocracy_bad <- wvs$Q239


#theocracy bad
wvs$expertrule_bad <- wvs$Q236

#index
wvs$democ_index <- rowSums(cbind(wvs$democracy_good, wvs$armyrule_bad, wvs$strongleader_bad, wvs$theocracy_bad, wvs$expertrule_bad))

sd(wvs$democ_index, na.rm=T)

###################################################################
####################        CONTROLS     ##########################
###################################################################

wvs$educ <- wvs$Q275

wvs$income_satisfied <- wvs$Q50

wvs$employed <- wvs$Q279

wvs$sex[wvs$Q260==2] <- 1
wvs$sex[wvs$Q260==1] <- 0


###ONLY WITH THE INDEX (TABLE 7.5)
summary(mod1 <- lm(democ_index ~ natharbdist_km + as.factor(B_COUNTRY)-1, data=wvs))
summary(mod2 <- lm(democ_index ~ natharbdist_km + urban + as.factor(B_COUNTRY)-1, data=wvs))
summary(mod3 <- lm(democ_index ~ natharbdist_km + educ + income_satisfied + sex + employed + urban + as.factor(B_COUNTRY)-1, data=wvs))

cluster_se1 <- as.vector(summary(mod1,cluster = c("gid"))$coefficients[,"Std. Error"])
cluster_se2 <- as.vector(summary(mod2,cluster = c("gid"))$coefficients[,"Std. Error"])
cluster_se3 <- as.vector(summary(mod3,cluster = c("gid"))$coefficients[,"Std. Error"])


library(stargazer)
# print stargazer output with robust standard errors
stargazer(mod1,mod2,mod3,
          type = "text",se = list(cluster_se1, cluster_se2,
                                  cluster_se3), out="output/table_7_5.txt")
